Ground-state of graphene in the presence of random charged impurities. 
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We calculate the carrier density dependent ground state properties of graphene in the presence 
of random charged impurities in the substrate taking into account disorder and interaction ef- 
fects non-perturbatively on an equal footing in a self-consistent theoretical formalism. We provide 
detailed quantitative results on the dependence of the disorder-induced spatially inhomogeneous 
two-dimensional carrier density distribution on the external gate bias, the impurity density, and the 
impurity location. We find that the interplay between disorder and interaction is strong, particu- 
larly at lower impurity densities. We show that for the currently available typical graphene samples, 
inhomogeneity dominates graphene physics at low (< 10 12 cm -2 ) carrier density with the density 
fluctuations becoming larger than the average density. 

PACS numbers: 



The recent experimental realization [jO] of single-layer 
graphene sheets has spurred an enormous amount of ac- 
tivity in studying the electronic properties of 2D chiral 
Dirac fermions in the context of solid state materials 
physics. While much of this interest is fundamental, a 
substantial part of it also derives from the technologi- 
cal prospect of graphene being used as a novel transis- 
tor material. To understand current experiments and 
be able to design future graphene-based electronic de- 
vices it is essential to know the properties, origin and 
effects of extrinsic disorder in graphene. The low energy 
electronic states of graphene are described by a massless 
Dirac equation. In clean isolated graphene (the so-called 
intrinsic graphene) the Fermi energy lies exactly at the 
Dirac point (i.e. the charge neutrality point) where the 
linear chiral electron and hole bands cross each other. 
Several works Q calculated the graphene conductivity 
assuming the graphene Fermi energy to be exactly at the 
Dirac point throughout the graphene layer. These works 
found the Dirac point conductivity to be either 0, oo or, 
in the limit of vanishing disorder, equal to the universal 
value <jd = 4e 2 /7rS. In current experiments however the 
measured conductivity at the Dirac point Q is finite and 
much bigger (by a factor of 2-20) than the universal the- 
oretical prediction <jd and varies strongly from sample to 
sample. The discrepancy can be resolved if we consider 
that disorder in addition to represent the main source of 
scattering has an other important effect: it locally shifts 
the Dirac point removing, at zero gate voltage, the Fermi 
energy from the charge neutrality point [4|. This leads 
immediately to a disorder-induced inhomogeneous den- 
sity landscape with electron-hole puddles. Such puddles 
have been proposed theoretically [5j and observed exper- 
imentally [6j, lZ|. Experiments, by themselves, are unable 
to directly identify the cause of the carrier density inho- 
mogeneities. Two kinds of disorder have been proposed 
in graphene to have this effect: ripples Q and random 
charge impurities 0. Transport theories [H, @, [l(J U| 
based on the presence of charge impurities have been 



successful in explaining the experimental results [3| . But 
whether the puddles arise from the random charged im- 
purities or from some other mechanism [8( has remained 
an open question. We provide in this Letter the first re- 
alistic theoretical description of the electron-hole puddles 
in graphene assuming the random charged impurity dis- 
order to be the underlying mechanism. Our theoretical 
results are in excellent qualitative agreement with the 
existing experimental data 0, 0|- A quantitative com- 
parison between our results and future experiments with 
higher quantitative accuracy would enable a definitive 
understanding of the nature of the disorder in graphene. 

At low energies the quasiparticles in graphene can be 
described by a massless Dirac- fermion, MDF, model with 
an ultraviolet cutoff wave vector k c . We set k c — 1/ao, 
where ao is the graphene lattice constant, ao = 0.246 nm, 
corresponding to an energy cut-off E c « 3 eV, and mea- 
sure the energies from the Dirac point. To find the 
ground-state carrier density n we use the Thomas-Fermi- 
Dirac (TF) theory. In contrast with the standard TF the- 
ory we retain the exchange potential non perturbatively 
through its local density approximation so that the en- 
ergy functional E[n] reads: 
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where Vf = 10 6 m/s is the Fermi velocity, r s = e 2 /(HvFe) 
is the coupling constant with e the effective background 
dielectric constant, ^cfn] is the exchange energy, Vd is 
the disorder potential and \i is the chemical potential. 
The first two terms in ([1]) are the kinetic energy and 
the Hartree part of the Coulomb interaction, respectively. 
For graphene on SiC>2 substrate e = 2.5 and then r s = 
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0.8. By differentiating E\n\ with respect to n we find: 
5E 



nxio'icm f < b > 
6000 r 



5n 



-flVf 



/ \ n— r »"s fn(r')cPr' 
sgn(n) ^ / ^ +r s V^ 



-S(n) — 
(2) 

where is the Hartree-Fock seif-energy [L2|,|l4| eval- 
uated at the Fermi wave vector kp = sgn(n)i/ 7r|n|: 
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where C « 0.916. 

We assume Vd to be the 2D Coulomb potential in the 
graphene plane generated by a random 2D distribution, 
C(r), of impurity charges placed at a distance ci from 
the graphene layer. Denoting by angular brackets the 
average over disorder realizations we assume: 



<C(r)> = 0; (C(n)C(r 2 )) = n imp <5(r 2 - n); (4) 



where 



L imp 



is the 2D impurity density. A non zero value 



of (C(r)) can be taken into account by a shift of fi, i.e. 
of the gate voltage. n imp and d should be taken as effec- 
tive parameters characterizing the impurity distribution 
in a minimal two parameter model. In current graphene 
samples obtained through mechanical exfoliation possible 
sources of charge impurities are most likely ions in the 
substrate that drift close to the surface, charges trapped 
between the graphene layer and the substrate and free 
charges that stick to the top surface of the graphene layer. 
This picture is consistent with the vast literature on dis- 
order in Si MOSFET and has recently been indir ectly 
confirmed by experiments on suspended graphene [15| . 
The values of n imp extracted from transport measure- 
ments, and used in this work, are indeed of the same order 
of magnitude, [10 11 — 10 12 ] cm -2 , as the ones used to de- 
scribe quantitatively disorder effects of MOSFET devices 
on SiC>2. Combining Eq.([2])-Q we find the ground state 
carrier density by solving the equation SE/dn — using 
the steepest descent method. Our calculations are done 
for a finite square lattice of size LxL. All the results pre- 
sented in this Letter are obtained for L = 200 nm and are 
found to be independent of system size for L > 100 nm. 
For the discretization in real space we use a 1 nm step. 

For a given disorder realization, for fi = 0, a typical re- 
sult, including exchange, for n(r) is shown in Fig. [1] (a). 
For n > (ra < 0) we have particles (holes). The result 
without exchange is characterized by larger density fluc- 
tuations. This is clear from Fig. [1] (b) which shows that 
the density distribution is more strongly peaked around 
n = when exchange is taken into account. The re- 
sult of Fig. Q] (b) is counter-intuitive because exchange 
suppresses density inhomogeneity instead of enhancing it 
as in parabolic-band inhomogencous electron liquids. A 
complementary DFT-LDA calculation, using single disor- 
der realizations with few impurities, has also found simi- 
lar results [l3[ . Contrary to our work in [l3| the correla- 
tion contributions have been taken into account. Given 
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FIG. 1: (Color online). Results at the Dirac point for a 
disorder realization assuming rii mp — 10 12 cm~ 2 , d = 1 nm, 
e = 2.5. (a) Color plot of n(r) including exchange, (b) Den- 
sity distribution for n(r) shown in (a). For clarity the result 
without exchange has been offset along the i-axis. 



the numerical complexity of the DFT-LDA approach in 
[l3T | only small samples were considered and disorder av- 
eraged results, that would permit a close quantitative 
comparison, were not presented. The results for sin- 
gle disorder realizations are qualitatively similar to ours, 
showing that correlation terms have only a minor quanti- 
tative effect. The reason is that in graphene to very good 
approximation the correlation term scales with n in the 
same way as exchange 12J, 1 1 31 ] but with opposite sign 
and therefore its effect is to simply reduce the exchange 
strength. 

The results of Fig. Q] are visually very similar to the 
ones observed in experiments @, 0], but a quantita- 
tive comparison can only be achieved by calculating 
the disorder-averaged statistical properties. For a given 
quantity X we therefore calculate its disorder averaged 
value, (X), and spatial correlation function ((SX(r)) 2 ) — 
{(X(r) - (X))(X(Q) - (X))). From these re sults we ex- 
tract the rms of the fluctuations, X rms = \J ((5X(0)) 2 ), 
and their typical correlation length £x = FWHM of 
((SX(r)) 2 ). At the neutrality point £ = £„ can be loosely 
taken as a measure of the electron-hole puddle size. In 
Fig. we present the disordered averaged results at the 
Dirac point as a function of ni mp . We see that exchange 
suppresses the amplitude of the density fluctuations and 
increases their correlation length and that its effect be- 
comes increasingly important as the impurity density de- 



creases; for the lowest ni mp the value of n rms includ- 
ing exchange is 3 times smaller than the value obtained 
without exchange, Fig. [2] (a). In addition we see that 
the scaling of n rms with nj mp is very different with and 
without exchange. From Fig. [2] (b) we see that as rii mp 
decreases £ increases very slowly, especially for low val- 
ues of d, a result that underlines the importance of non- 
linear screening terms. Adapted to a 2D distribution of 
charges the approach used in Ref. [l6| for the scaling of 
£ on n lmp gives £ s» l/{r 2 ^/n lmp ). For r s = 0.8 and 
nimp = 2. 10 9 cm -2 we would then expect £ m 350 nm, a 
value an order of magnitude larger than the value shown 
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in Fig. [2] (b). The reason for this discrepancy is that for 
small values of rii mp the carrier distribution is not char- 
acterized by smooth long-range fluctuations but rather 
by wide regions of very small carrier density (« 0) inter- 
spersed with small electron-hole puddles with the typical 
size £ shown in Fig. [2] (b). This picture is confirmed in 
Fig. [2(c) where the disorder averaged area fraction, Aq, 
over which |n(r) — (n)| < n rms /10 is plotted as a func- 
tion of rii m p. We see that as rn mp decreases A increases 
reaching more than 1/3 at the lowest impurity densities. 
The fraction of area over which |n(r) — (n)| is less than 
1/5 of n rms surpasses 50% for rii mp < 10 10 cm -2 . Thus, 
much of the 2D landscape in this situation has very low 
(<C rii mp ) carrier density with a few random electron-hole 
puddles. In Fig.[2](d) the dependence of the excess charge 
SQ = n rms Tr^ 2 on n,i mp is shown. At high impurity den- 
sities (> 10 12 cm" 2 ) and values of d > 1 nm, SQ can be 
approximately identified with the average number of car- 
riers per puddle, however the above discussion and the re- 
sults for Aq allow us to recognize that SQ, for small TH mp , 
is not the typical number of carriers per puddle. The 
reason is that for small rii mp (and/or d), because of the 
large fraction of area over which is \n(r) — (n)| <C n rms , 
n rm s is much smaller than the typical carrier density in 
a electron-hole puddle of size £. At low rii mp SQ grossly 
underestimates the number of carriers in a typical puddle 
of size £. 




FIG. 2: (Color online). Results as a function of rn mp at the 
Dirac point for d — 1 nm (blue lines) and d = 0.3 nm (red 
lines), e = 2.5. (a) n rms ; (b) £ ; (c) A ; (d) SQ. 



We are now in a position to discuss the validity of our 
TF approach. The use of the TF theory is justified when 
the condition \\7n(r)\/[n(r)kp(r)) <C 1 is satisfied. If we 
estimate |Vn(r)| w n rms /^, the above inequality implies 
y/irn rrns t; ^> 1, i.e. SQ 3> 1. However, for small ni mp and 
d, n rrns greatly underestimates |n| in the regions where 
it is inhomogeneous, i.e. in the electron-hole puddles of 



size £. We find that at low n,i mp , \n\ in these electron- hole 
puddles is a factor of 10 or more higher than n rms . This 
can already be seen for relatively high values of ni mp and 
d: from Fig.Q](a) we see that \n\ inside the electron-hole 
puddles takes values as high as 8. 10 12 cm" 2 whereas 
the corresponding value of n rms is only 8. 10 11 cm" 2 , 
Fig. [2] (a) . Even in the limiting case of an isolated impu- 
rity with d = the density profile obtained using the TF 
approach [13] is very similar to the one obtained starting 
from the Dirac equation and treating the interaction via 
RG [l|. The additional 8(r) for n(r) found in (and 
[l^l ) in real graphene, in which the MDF model ap plie s 
only at low energies, is regularized by max[d,ao], [201 ] . 
Our results are therefore quantitatively accurate. 

From the theoretical analysis 0, 0, 0, EH of ex- 
perimental transport results [§] one obtains for typi- 
cal graphene samples on SiC>2, d = 1 nm and n s» 
3. 10 11 cm " 2 . For these values, from Fig. [2] (a) and 
(b) we see that n rms — 3. 10 11 cm ~ 2 and £ = 9 nm. The 
value of n rms is in very good agreement with the recent 
STM and SET @ results. The value of £ is also in 
very good agreement with the STM results and consis- 
tent with the results of Ref. Q that, given the lower SET 
spatial resolution (> 150 nm), could only provide for £ 
an upper bound of 30 nm. 

At a finite gate voltage, V g , the average carrier density 
(n) = CgVg/e is induced, where C g is the gate capaci- 
tance. In our calculations we indirectly fix (n) by vary- 
ing the chemical potential /i. The relation between fj, and 
(n) is shown in Fig. [3] (a). Contrary to ordinary parabolic 
band fermionic systems the relation between fj, and (n) 
strongly depends on disorder even when exchange is ne- 
glected. This is also shown in Fig. [3] (b) in which /i is 
plotted as a function of rii mp for a fixed value of (n) . The 
dependence of fJ,((n)) on rii mp even without exchange is 
due to the fact that in graphene the kinetic energy does 
not scale linearly with n. From Fig. [3] (a) we see that only 
for riimp < 10 10 cm" 2 /i((n)) follows the equation valid 
for clean graphene. We also notice that n((n)) is strongly 
affected by exchange. The results of Fig. [3] demonstrate 
the interplay of disorder and interaction in graphene and 
show how the dependence of /1 on (n), and in particular 
the average compressibility l/(n 2 dn/dn), can be used to 
probe the strength of disorder and many-body effects. 
In Fig. [5] (a) the disorder averaged density distribution 
obtained including exchange is plotted for different val- 
ues of (n), i.e. of V g . For V g — the distribution has 
a strong peak (ss 20 times the maximum of the y-scale 
of Fig. H] (a)) at n — 0. As V g increases the n = 
peak survives and a broad peak at finite n develops. For 
large enough V g the n = peak disappears and the den- 
sity distribution is characterized only by the broad peak 
centered at n = (ri). The results without exchange are 
qualitatively similar. The double peak structure for finite 
V g provides direct evidence for the existence of puddles 
over a finite voltage range. High values of V g remove one 
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FIG. 3: (Color online). Solid (dashed) lines: results with 
(without) exchange, e = 2.5. (a) fi vs. (n) for rii mp = 
2.10 10 cm -2 , red lines; rii mp — 10 12 cm -2 , blue lines; 
riimp = 2.10 12 cm -2 , green lines. The black dot-dash (dash) 
line shows n(n) for clean graphene with (without) exchange, 
(b) n vs. n imp for (n) = 5 10 11 cm -2 . 



small puddles and large regions of almost zero (|n| <§; 
nimp) carrier density; (4) in current samples n rms > (n) 
for carrier densities as high as (n) 10 12 cm~ 2 ; (5) the 
number of carriers per puddle is ^1-5 at low carrier den- 
sities; (6) our theory agrees well with the existing data 
0, 0] but more experiments will be required to test our 
quantitative predictions. 
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kind of puddles and increase the amplitude of the den- 
sity fluctuations reflected in an increase of n rms . On the 
other hand the ratio n rms /(n) decreases monotonically 
as a function of (n) as can be seen in Fig. [4] (b). We 
can define a characteristic density, the value of (n) 

for which n rms — (n) with AV g — en c /C g loosely mea- 
suring the width in gate voltage over which the transport 
properties of graphene are dominated by the density fluc- 
tuations around the Dirac point. The inset of Fig. [4] (b) 
shows n c as a function of n; m „ for d = 1 nm. We can see 
that in current samples n 
as high as (n) « 10 
of the carrier density distribution and n rms on V g are 
unique to inhomogeneities created by charged impurities 
and is a prediction that should be easy to verify experi- 
mentally 
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FIG. 4: (Color online). Results away from the Dirac point 
assuming Si02 substrate, (a) Density distribution averaged 
over disorder for different values of V g for d = 1 nm and 
riimn = 10 12 cm~ 2 . (b) n rms /(n) vs. (n) for d = 1 nm and: 
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= 1.5 x 10 cm , red lines; 



n 



imp 



blue 



lines; 



5 x 10 cm , green lines. Inset: n c vs. 



for d = 1 nm. Solid (dashed) lines: results with (without) 
exchange. 



We conclude by summarizing our key qualitative find- 
ings: (1) both disorder and many-body effects become 
quantitatively very important on the chemical potential 
close to the Dirac point; (2) many body effects are more 
important at lower values of ni mp ; (3) for low rii mv the 
ground state near the Dirac point is characterized by 
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